﻿/*
Adapted from AAPlus library by PJ Naughter
Permission was granted to convert the library to C#

This source is modified from PJ Naugters AA+ lib
(http://www.naughter.com/aa.html)
and therefore is covered under his copyright license, which
is as follows:

Purpose: Implementation for the algorithms which obtain the physical parameters of the Moon
Created: PJN / 17-01-2004
History: PJN / 19-02-2004 1. The optical libration in longitude is now returned in the range -180 - 180 degrees

Copyright (c) 2004 - 2009 by PJ Naughter (Web: www.naughter.com, Email: pjna@naughter.com)

All rights reserved.

Copyright / Usage Details:

You are allowed to include the source code in any product (commercial, shareware, freeware or otherwise) 
when your product is released in binary form. You are allowed to modify the source code in any way you want 
except you cannot modify the copyright details at the top of each module. If you want to distribute source 
code with your application, then you are only allowed to distribute versions released by the author. This is 
to maintain a single distribution point for the source code. 

*/


using System;

namespace AAPlus
{
    public class PhysicalMoonDetails
    {
        public double ldash;
        public double bdash;
        public double ldash2;
        public double bdash2;
        public double l;
        public double b;
        public double P;
    }

    public class SelenographicMoonDetails
    {
        public double l0;
        public double b0;
        public double c0;

        public SelenographicMoonDetails()
        {
        }

        public SelenographicMoonDetails(double fl0, double fb0, double fc0)
        {
            l0 = fl0;
            b0 = fb0;
            c0 = fc0;
        }
    };

    public static class PhysicalMoon
    {


        public static void CalculateOpticalLibration(double JD, double Lambda, double Beta, out double ldash, out double bdash, out double ldash2, out double bdash2, out double epsilon, out double omega, out double DeltaU, out double sigma, out double I, out double rho)
        {
            //Calculate the initial quantities
            double Lambdarad = CoordinateTransformation.DegreesToRadians(Lambda);
            double Betarad = CoordinateTransformation.DegreesToRadians(Beta);
            I = CoordinateTransformation.DegreesToRadians(1.54242);
            DeltaU = CoordinateTransformation.DegreesToRadians(Nutation.NutationInLongitude(JD) / 3600);
            double F = CoordinateTransformation.DegreesToRadians(Moon.ArgumentOfLatitude(JD));
            omega = CoordinateTransformation.DegreesToRadians(Moon.MeanLongitudeAscendingNode(JD));
            epsilon = Nutation.MeanObliquityOfEcliptic(JD) + Nutation.NutationInObliquity(JD) / 3600;

            //Calculate the optical librations
            double W = Lambdarad - DeltaU / 3600 - omega;
            double A = Math.Atan2(Math.Sin(W) * Math.Cos(Betarad) * Math.Cos(I) - Math.Sin(Betarad) * Math.Sin(I), Math.Cos(W) * Math.Cos(Betarad));
            ldash = CoordinateTransformation.MapTo0To360Range(CoordinateTransformation.RadiansToDegrees(A) - CoordinateTransformation.RadiansToDegrees(F));
            if (ldash > 180)
                ldash -= 360;
            bdash = Math.Asin(-Math.Sin(W) * Math.Cos(Betarad) * Math.Sin(I) - Math.Sin(Betarad) * Math.Cos(I));

            //Calculate the physical librations
            double T = (JD - 2451545.0) / 36525;
            double K1 = 119.75 + 131.849 * T;
            K1 = CoordinateTransformation.DegreesToRadians(K1);
            double K2 = 72.56 + 20.186 * T;
            K2 = CoordinateTransformation.DegreesToRadians(K2);

            double M = Earth.SunMeanAnomaly(JD);
            M = CoordinateTransformation.DegreesToRadians(M);
            double Mdash = Moon.MeanAnomaly(JD);
            Mdash = CoordinateTransformation.DegreesToRadians(Mdash);
            double D = Moon.MeanElongation(JD);
            D = CoordinateTransformation.DegreesToRadians(D);
            double E = Earth.Eccentricity(JD);

            rho = -0.02752 * Math.Cos(Mdash) +
                  -0.02245 * Math.Sin(F) +
                  0.00684 * Math.Cos(Mdash - 2 * F) +
                  -0.00293 * Math.Cos(2 * F) +
                  -0.00085 * Math.Cos(2 * F - 2 * D) +
                  -0.00054 * Math.Cos(Mdash - 2 * D) +
                  -0.00020 * Math.Sin(Mdash + F) +
                  -0.00020 * Math.Cos(Mdash + 2 * F) +
                  -0.00020 * Math.Cos(Mdash - F) +
                  0.00014 * Math.Cos(Mdash + 2 * F - 2 * D);

            sigma = -0.02816 * Math.Sin(Mdash) +
                    0.02244 * Math.Cos(F) +
                    -0.00682 * Math.Sin(Mdash - 2 * F) +
                    -0.00279 * Math.Sin(2 * F) +
                    -0.00083 * Math.Sin(2 * F - 2 * D) +
                    0.00069 * Math.Sin(Mdash - 2 * D) +
                    0.00040 * Math.Cos(Mdash + F) +
                    -0.00025 * Math.Sin(2 * Mdash) +
                    -0.00023 * Math.Sin(Mdash + 2 * F) +
                    0.00020 * Math.Cos(Mdash - F) +
                    0.00019 * Math.Sin(Mdash - F) +
                    0.00013 * Math.Sin(Mdash + 2 * F - 2 * D) +
                    -0.00010 * Math.Cos(Mdash - 3 * F);

            double tau = 0.02520 * E * Math.Sin(M) +
                         0.00473 * Math.Sin(2 * Mdash - 2 * F) +
                         -0.00467 * Math.Sin(Mdash) +
                         0.00396 * Math.Sin(K1) +
                         0.00276 * Math.Sin(2 * Mdash - 2 * D) +
                         0.00196 * Math.Sin(omega) +
                         -0.00183 * Math.Cos(Mdash - F) +
                         0.00115 * Math.Sin(Mdash - 2 * D) +
                         -0.00096 * Math.Sin(Mdash - D) +
                         0.00046 * Math.Sin(2 * F - 2 * D) +
                         -0.00039 * Math.Sin(Mdash - F) +
                         -0.00032 * Math.Sin(Mdash - M - D) +
                         0.00027 * Math.Sin(2 * Mdash - M - 2 * D) +
                         0.00023 * Math.Sin(K2) +
                         -0.00014 * Math.Sin(2 * D) +
                         0.00014 * Math.Cos(2 * Mdash - 2 * F) +
                         -0.00012 * Math.Sin(Mdash - 2 * F) +
                         -0.00012 * Math.Sin(2 * Mdash) +
                         0.00011 * Math.Sin(2 * Mdash - 2 * M - 2 * D);

            ldash2 = -tau + (rho * Math.Cos(A) + sigma * Math.Sin(A)) * Math.Tan(bdash);
            bdash = CoordinateTransformation.RadiansToDegrees(bdash);
            bdash2 = sigma * Math.Cos(A) - rho * Math.Sin(A);
        }

        public static PhysicalMoonDetails CalculateHelper(double JD, out double Lambda, out double Beta, out double epsilon, out Coordinate2D Equatorial)
        {
            //What will be the return value
            PhysicalMoonDetails details = new PhysicalMoonDetails();

            //Calculate the initial quantities
            Lambda = Moon.EclipticLongitude(JD);
            Beta = Moon.EclipticLatitude(JD);

            //Calculate the optical libration
            double omega;
            double DeltaU;
            double sigma;
            double I;
            double rho;
            CalculateOpticalLibration(JD, Lambda, Beta, out details.ldash, out details.bdash, out details.ldash2, out details.bdash2, out epsilon, out omega, out DeltaU, out sigma, out I, out rho);
            double epsilonrad = CoordinateTransformation.DegreesToRadians(epsilon);

            //Calculate the total libration
            details.l = details.ldash + details.ldash2;
            details.b = details.bdash + details.bdash2;
            double b = CoordinateTransformation.DegreesToRadians(details.b);

            //Calculate the position angle
            double V = omega + DeltaU + CoordinateTransformation.DegreesToRadians(sigma) / Math.Sin(I);
            double I_rho = I + CoordinateTransformation.DegreesToRadians(rho);
            double X = Math.Sin(I_rho) * Math.Sin(V);
            double Y = Math.Sin(I_rho) * Math.Cos(V) * Math.Cos(epsilonrad) - Math.Cos(I_rho) * Math.Sin(epsilonrad);
            double w = Math.Atan2(X, Y);

            Equatorial = CoordinateTransformation.Ecliptic2Equatorial(Lambda, Beta, epsilon);
            double Alpha = CoordinateTransformation.HoursToRadians(Equatorial.X);

            details.P = CoordinateTransformation.RadiansToDegrees(Math.Asin(Math.Sqrt(X * X + Y * Y) * Math.Cos(Alpha - w) / (Math.Cos(b))));

            return details;
        }

        public static PhysicalMoonDetails CalculateGeocentric(double JD)
        {
            double Lambda;
            double Beta;
            double epsilon;
            Coordinate2D Equatorial = new Coordinate2D();
            return CalculateHelper(JD, out Lambda, out Beta, out epsilon, out Equatorial);
        }

        public static PhysicalMoonDetails CalculateTopocentric(double JD, double Longitude, double Latitude)
        {
            //First convert to radians
            Longitude = CoordinateTransformation.DegreesToRadians(Longitude);
            Latitude = CoordinateTransformation.DegreesToRadians(Latitude);

            double Lambda;
            double Beta;
            double epsilon;
            Coordinate2D Equatorial = new Coordinate2D();
            PhysicalMoonDetails details = CalculateHelper(JD, out Lambda, out Beta, out epsilon, out Equatorial);

            double R = Moon.RadiusVector(JD);
            double pi = Moon.RadiusVectorToHorizontalParallax(R);
            double Alpha = CoordinateTransformation.HoursToRadians(Equatorial.X);
            double Delta = CoordinateTransformation.DegreesToRadians(Equatorial.Y);

            double AST = Sidereal.ApparentGreenwichSiderealTime(JD);
            double H = CoordinateTransformation.HoursToRadians(AST) - Longitude - Alpha;

            double Q = Math.Atan2(Math.Cos(Latitude) * Math.Sin(H), Math.Cos(Delta) * Math.Sin(Latitude) - Math.Sin(Delta) * Math.Cos(Latitude) * Math.Cos(H));
            double Z = Math.Acos(Math.Sin(Delta) * Math.Sin(Latitude) + Math.Cos(Delta) * Math.Cos(Latitude) * Math.Cos(H));
            double pidash = pi * (Math.Sin(Z) + 0.0084 * Math.Sin(2 * Z));

            double Prad = CoordinateTransformation.DegreesToRadians(details.P);

            double DeltaL = -pidash * Math.Sin(Q - Prad) / Math.Cos(CoordinateTransformation.DegreesToRadians(details.b));
            details.l += DeltaL;
            double DeltaB = pidash * Math.Cos(Q - Prad);
            details.b += DeltaB;
            details.P += DeltaL * Math.Sin(CoordinateTransformation.DegreesToRadians(details.b)) - pidash * Math.Sin(Q) * Math.Tan(Delta);

            return details;
        }

        public static SelenographicMoonDetails CalculateSelenographicPositionOfSun(double JD)
        {
            double R = Earth.RadiusVector(JD) * 149597970;
            double Delta = Moon.RadiusVector(JD);
            double lambda0 = Sun.ApparentEclipticLongitude(JD);
            double lambda = Moon.EclipticLongitude(JD);
            double beta = Moon.EclipticLatitude(JD);

            double lambdah = CoordinateTransformation.MapTo0To360Range(lambda0 + 180 + Delta / R * 57.296 * Math.Cos(CoordinateTransformation.DegreesToRadians(beta)) * Math.Sin(CoordinateTransformation.DegreesToRadians(lambda0 - lambda)));
            double betah = Delta / R * beta;

            //What will be the return value
            SelenographicMoonDetails details = new SelenographicMoonDetails();

            //Calculate the optical libration
            double omega;
            double DeltaU;
            double sigma;
            double I;
            double rho;
            double ldash0;
            double bdash0;
            double ldash20;
            double bdash20;
            double epsilon;
            CalculateOpticalLibration(JD, lambdah, betah, out ldash0, out bdash0, out ldash20, out bdash20, out epsilon, out omega, out DeltaU, out sigma, out I, out rho);

            details.l0 = ldash0 + ldash20;
            details.b0 = bdash0 + bdash20;
            details.c0 = CoordinateTransformation.MapTo0To360Range(450 - details.l0);
            return details;
        }

        public static double AltitudeOfSun(double JD, double Longitude, double Latitude)
        {
            //Calculate the selenographic details
            SelenographicMoonDetails selenographicDetails = CalculateSelenographicPositionOfSun(JD);

            //convert to radians
            Latitude = CoordinateTransformation.DegreesToRadians(Latitude);
            Longitude = CoordinateTransformation.DegreesToRadians(Longitude);
            selenographicDetails.b0 = CoordinateTransformation.DegreesToRadians(selenographicDetails.b0);
            selenographicDetails.c0 = CoordinateTransformation.DegreesToRadians(selenographicDetails.c0);

            return CoordinateTransformation.RadiansToDegrees(Math.Asin(Math.Sin(selenographicDetails.b0) * Math.Sin(Latitude) + Math.Cos(selenographicDetails.b0) * Math.Cos(Latitude) * Math.Sin(selenographicDetails.c0 + Longitude)));
        }

        public static double SunriseSunsetHelper(double JD, double Longitude, double Latitude, bool bSunrise)
        {
            double JDResult = JD;
            double Latituderad = CoordinateTransformation.DegreesToRadians(Latitude);
            double h;
            do
            {
                h = AltitudeOfSun(JDResult, Longitude, Latitude);
                double DeltaJD = h / (12.19075 * Math.Cos(Latituderad));
                if (bSunrise)
                    JDResult -= DeltaJD;
                else
                    JDResult += DeltaJD;
            }
            while (Math.Abs(h) > 0.001);

            return JDResult;
        }

        public static double TimeOfSunrise(double JD, double Longitude, double Latitude)
        {
            return SunriseSunsetHelper(JD, Longitude, Latitude, true);
        }

        public static double TimeOfSunset(double JD, double Longitude, double Latitude)
        {
            return SunriseSunsetHelper(JD, Longitude, Latitude, false);
        }
    }
}
